#include "sound-wave.hh"
#include <map>
#include <math.h>

namespace Dune {
  namespace Equation {
    SoundWave_<1>::SoundWave_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    SoundWave_<1>::~SoundWave_(){
    }
    SoundWave_<2>::SoundWave_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    SoundWave_<2>::~SoundWave_(){
    }
    SoundWave_<3>::SoundWave_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    SoundWave_<3>::~SoundWave_(){
    }

    void SoundWave_<1>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("velocity1");
      set_equ_para("u0",0.0);
      set_equ_para("rho0",0.0);
      set_equ_para("k0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["rho0"]=ptree.get("equation.rho0",0.0);
      para["k0"]=ptree.get("equation.k0",0.0);
    }
    void SoundWave_<2>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("velocity1");
      varname.push_back("velocity2");
      set_equ_para("u0",0.0);
      set_equ_para("v0",0.0);
      set_equ_para("rho0",0.0);
      set_equ_para("k0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["v0"]=ptree.get("equation.v0",0.0);
      para["rho0"]=ptree.get("equation.rho0",0.0);
      para["k0"]=ptree.get("equation.k0",0.0);
    }
    void SoundWave_<3>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("velocity1");
      varname.push_back("velocity2");
      varname.push_back("velocity3");
      set_equ_para("u0",0.0);
      set_equ_para("v0",0.0);
      set_equ_para("w0",0.0);
      set_equ_para("rho0",0.0);
      set_equ_para("k0",0.0);
      para["u0"]=ptree.get("equation.u0",0.0);
      para["v0"]=ptree.get("equation.v0",0.0);
      para["w0"]=ptree.get("equation.w0",0.0);
      para["rho0"]=ptree.get("equation.rho0",0.0);
      para["k0"]=ptree.get("equation.k0",0.0);
      // DEBUG
      std::cout << "equation init: u0="
        << para["u0"] << std::endl;
      std::cout << "equation init: v0="
        << para["v0"] << std::endl;
      std::cout << "equation init: w0="
        << para["w0"] << std::endl;
      std::cout << "equation init: rho0="
        << para["rho0"] << std::endl;
      std::cout << "equation init: k0="
        << para["k0"] << std::endl;
    }

    // FIXME: need improvement
    double SoundWave_<1>::local_max_eigenvalue(double* var){
      return  para["k0"]/para["rho0"]+para["u0"];
    }
    double SoundWave_<2>::local_max_eigenvalue(double* var){
      return para["k0"]/para["rho0"]+std::sqrt(std::pow(para["u0"],2)+std::pow(para["v0"],2));
    }
    double SoundWave_<3>::local_max_eigenvalue(double* var){
      return para["k0"]/para["rho0"]+std::sqrt(std::pow(para["u0"],2)+std::pow(para["v0"],2)+std::pow(para["w0"],2));   
    }

    void SoundWave_<1>::get_equation_flux(int whichdim, double* fequ, double* var){
      get_1d_fx(fequ,var);
    }
    void SoundWave_<2>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_2d_fx(fequ,var);
      if(whichdim==1)
        get_2d_fy(fequ,var);
    }
    void SoundWave_<3>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_3d_fx(fequ,var);
      if(whichdim==1)
        get_3d_fy(fequ,var);
      if(whichdim==2)
        get_3d_fz(fequ,var);
    }
    
    void SoundWave_<1>::get_1d_fx(double* fequ, double* var){
      fequ[0] = para["u0"]      *var[0]+para["k0"]*var[1];
      fequ[1] = 1.0/para["rho0"]*var[0]+para["u0"]*var[1];  
    }
    void SoundWave_<2>::get_2d_fx(double* fequ, double* var){
      fequ[0] = para["u0"]      *var[0]+para["k0"]*var[1];
      fequ[1] = 1.0/para["rho0"]*var[0]+para["u0"]*var[1];
      fequ[2] =                         para["u0"]*var[2];      
    }
    void SoundWave_<2>::get_2d_fy(double* fequ, double* var){
      fequ[0] = para["v0"]      *var[0]+para["k0"]*var[2];
      fequ[1] =                         para["v0"]*var[1];
      fequ[2] = 1.0/para["rho0"]*var[0]+para["v0"]*var[2];  
    }
    void SoundWave_<3>::get_3d_fx(double* fequ, double* var){
      fequ[0] = para["u0"]      *var[0]+para["k0"]*var[1];
      fequ[1] = 1.0/para["rho0"]*var[0]+para["u0"]*var[1];
      fequ[2] =                         para["u0"]*var[2];
      fequ[3] =                         para["u0"]*var[3];
    }
    void SoundWave_<3>::get_3d_fy(double* fequ, double* var){
      fequ[0] = para["v0"]      *var[0]+para["k0"]*var[2];
      fequ[1] =                         para["v0"]*var[1];
      fequ[2] = 1.0/para["rho0"]*var[0]+para["v0"]*var[2];
      fequ[3] =                         para["v0"]*var[3];    
    }
    void SoundWave_<3>::get_3d_fz(double* fequ, double* var){
      fequ[0] = para["w0"]      *var[0]+para["k0"]*var[3];
      fequ[1] =                         para["w0"]*var[1];
      fequ[2] =                         para["w0"]*var[2];
      fequ[3] = 1.0/para["rho0"]*var[0]+para["w0"]*var[3];
    }
  } // Equation
} // Dune
